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Global error estimation of linear 
multistep methods through the 
Runge-Kutta methods 


J. Farzi* 


Abstract 

In this paper, we study the global truncation error of the linear multi- 
step methods (LMM) in terms of local truncation error of the corresponding 
Runge-Kutta schemes. The key idea is the representation of LMM with a cor- 
responding Runge-Kutta method. For this, we need to consider the multiple 
step of a linear multistep method as a single step in the corresponding Runge- 
Kutta method. Therefore, the global error estimation of a LMM through the 
Runge-Kutta method will be provided. In this estimation, we do not take 
into account the effects of roundoff errors. The numerical illustrations show 
the accuracy and efficiency of the given estimation. 


Keywords: Linear multistep methods; Runge-Kutta methods; Local trun- 
cation error; Global error; Error estimation. 


1 Introduction 


The error estimation is one of the major issues in designing numerical algo- 
rithms. In the study of the linear multistep methods (LMM) 


k k 
Po Gling = Sine (1) 
jJ=0 j=0 
for solving an ordinary differential system 


ee a 


y(xo) = Yo; 


(2) 


where, f : Rx R™” — R”, there is a challenging issue of estimation of global 
truncation error (GTE) or simply global error. In spite the local trunca- 
tion error (LTE), the estimation of GTE is much more complicated. The 
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LTE estimations have been studied in some special cases, e.g., predictor- 
corrector (PC) and embedded Runge-Kutta methods. In the case of PC 
methods that predictor is an Adams-Bashforth scheme of order p and cor- 
rector is an Adams-Moulton scheme of the same order, the Milne estimation 
provides an estimation for LTE of the resulted PC method [9,10]. Recently, 
Cao and Petzold have developed an estimation for global error with adjoint 
method [3]. However, many efforts have previously been reported for provid- 
ing GTE bounds [6,7]. These bounds are of no practical value and the LTE 
form estimation of GTE, given in this paper, is more simpler than the usual 
theoretical bounds. For more extensive discussion on linear multistep meth- 
ods and their important subclauses, like backward difference formula (BDF) 
schemes see [1, 2,8-10, 12]. The more important application of the LMMs 
is in the time discretization of time dependent partial differential equations. 
The LMMs with strong stability preserving property (SSP) have a major role 
in this context [4,5]. 


In this paper, we represent a multiple step of a LMM as a single step of a 
new Runge-Kutta method. Then, we accomplish the GTE estimation of the 
LMM by estimation of LTE of the corresponding new Runge-Kutta method. 


This paper has been organized as follow: In Section 2 a very short review 
of the Runge-Kutta schemes is presented. Then, in Section 3 the main idea 
of the paper is given for one-step methods with a rather detailed study of the 
LTE and stability region of the new method. In Section 4 the idea of previous 
section is generalized to formulate a general LMM in the form of a new 
Runge-Kutta method. The new RK method has a popular structure in view 
of Butcher array. We take into account a starting procedure, subsequently, 
the method order and its LTE determined, which is the GTE of the original 
scheme. Finally, in Section 5 we present some numerical tests including a 
fourth order total variation bounded (TVB) scheme to illustrate efficiency of 
the given theory. 


2 A review of Runge-Kutta methods 


In this section, we shortly review the main concepts of Runge-Kutta methods 
that are required in the rest of the paper. For more details we refer to [2,9]. 
The s-stage Runge-Kutta method for the problem (2) is defined as follow 


Unt =Yn th S_ biki, (3) 
i=l 
where 


ki =f(tn+cihyyn +h > ashy), %=1,2,...,8. (4) 


j=l 
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Traditionally the Runge-Kutta methods is represented by the following 
Butcher array 


c|A 
o (6) 
where, 
C= Citra. sta]? ¢ B= Dig Basg bal A= lag a: (6) 


The derivation of a Runge-Kutta method of an arbitrary order is a crucial 
work without using the advanced concepts of elementary differentials and 
most related rooted trees. In fact, there is a 1-1 corresponding between ele- 
mentary differential (that are defined by Fréchet derivative) and the rooted 
trees. Therefore, using the rooted trees one can easily define the order con- 
dition and LTE of a Runge-Kutta method. The local truncation error of the 
p' order Runge-Kutta method (5) is given by 


Lre= YS ay -rvolro+owr), 
* r(t)=p+1 
where, () 
r(t)! 
Or) 


and r(t),a(t) and y(t) are order, symmetry and density of a tree t. The 
function F'(t) is defined on the set T of all trees which corresponds between 
the rooted trees and elementary differentials. The elementary differentials 
are evaluated at the value y(x,,). The w(t) function is also defined on the set 
of all trees T. For example we have, 


rh = =7, 


EM) a hig eeialtina: 


F(L) = {of}e = FY FOF), 


where, f()(K1, K2,..., Ky), K, € R™,t = 1,2,...,M is the Mth order 
Fréchet derivative of f. For more detailed definition of these functions see 
(2, 9]. 
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The following theorem gives the coefficients of the linear combination of 
y for a general q in terms of elementary differentials [9]: 


Theorem 1. Let y be the solution of the autonomous problem (1). Then 


y% = So a(t)F). (8) 
r(t)=q 
The stability function of a Runge-Kutta method is given by 


det [I — hA + heb™] 
det [I — hA] 


R(h) = (9) 


where, h= Ah, here is typically an eigenvalue of the jacobian matrix of f 
that is equivalently the eigenvalues of the linearized equation. 


3 One step methods 


In this section, we firstly study the situation for the one-step methods. The 
idea will be extended to a general linear multistep methods in the next sec- 
tions. Consider the following linear one step method 


Yn41 = Yn + h(Bofn + B1fn+1) (10) 


for a consistent method we have $9 + 6, = 1. Applying the above rule on N 
successive steps to advance the solution from x9 to xy we obtain 


yi = yo + h(Bofo + Pr fi) 
Yo = 41 + h(Bofi + Fife) 


yn = yn-1 + h(Bofn-1 + Af). 
Introducing the following slopes 


ky = f (0, Yo), ke = Piya kn = f(an,yn), 


we can write (10) as a (N+1)-stage Runge-Kutta method with the steplength 
H=Nh, 


HH 
Ym+1 = Ym + vy (bok: +ko+--:+kn + Bikn41) 


where, tm = Xo,tm41 = tm +H = xy and, 
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ky = f (tins Ym) 
ko = f(tmt an + (ok + Bi k2)) 


; H 
ky4i = f(tm + By ym + y (Bok +ko+---+kyw + Bikn41)). 


Thus, the corresponding Butcher array takes the following form 


0] 0 
1} Bo Br 

N|N N 
2/fo 1 fr 
N|N N N 

(11) 
N|/6o 1 .,, 1 fi 
N N N N 

Bo 1 A ft 

N N N 


3.1 Local truncation error 


In this section we obtain the LTE and the order of the reduced method 
(11). To determine the order of the deduced RK scheme we verify the order 
conditions for the Butcher array (11): 


N+1 
o(e)= 0 =1, (12) 
wa 


wd) = De bici = 5 (13) 


Since we have 6p + 3; = 1, substituting the data from (11) we observe 
that the first condition always holds 


N?-N+26.N_1 1 
= 2N2 “Sg EE: 
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where, we have used {9 + 8; = 1. Therefore, the only possible second order 
one step method is the trapezoidal rule. 


Now, we can find the general form of the local truncation error for both 
first (forward and backward euler) and second order (Trapezoidal) methods. 


For the second order method we have, 8; = $, 


N41 
v( eS ia Do eee 
i=1 


6N2 
N+1 2 
Pal 1 2N2 41 
t2=4) = D> biayey = <[z,1,...,1, <]4e= 
V(t =s) 2 oie ws gl4e = ona 


It is evident that in the limit the above order conditions, as N — ov, tend 
to the following infinite dimensional exact order conditions: 


The corresponding elementary differentials with the rooted trees ¢; and 
to are, 


F(ti) =f Ff), F(t) = fP F(A). 


To specify the LTE of second order scheme we need a representation of y‘?) 
in terms of elementary differentials. According to Theorem 1, we have 


yO = F(t,) + F(te), (14) 


therefore, the principle term in local truncation error (PLTE) for 6, = s, 
where p = 2 reads 


3 
PLIE = = Yoel — OVE (15) 
r(t)=3 
He 1 
= F(-aqa) (Ft) + Fe) 
= -~SNhiy® (ao) 
= - 5 len — 20)y (eo), 


and then, 
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LTE = —h" (an — £o)y (#o) + O(h3). (16) 


We can also regard (16) as the global error in N steps of the trapezoidal rule. 
Similarly, for 6, # 4 we obtain 


LTE = SN(N — 14 261)h?y (a) + O(A?). (17) 


3.2 Stability regions 


To construct the stability function of (11) we simply note that 


Bol...1 A 
i H | Go fi x Al ol 1h 
HA=— Heb? = — 
N eee al) 
Bo Te snsa 1.” Be Bol...1 AB 
thus we have 
‘7 Bo A 
1+Hy Nw N WPA 
1+H% H #8, 
I-HA+Heb™ = Be : : 
1 
and thereby, 
det{I — HA+ Heb™] = (1+ aN 


Similarly, we can show that the following relation is also valid 
say as eer 
det(I — HA] = (1- H—)”. 
N 
Inserting the above results into (9) we obtain 


wy (L+H) 


R(H) = a= Fay (18) 


which is the stability function of (11). 
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Stability region for By=1, B,=0 


4 T 


scaled by N 
oO 


-2 
scaled by N 
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Figure 1: The stability region of the new RK method with 89 = 1, 8; = 0 


Stability region for By=0, B,=1 


scaled by N 


Figure 2: The stability region of the new RK method with 69 = 0,6; =1 


1 2 
scaled by N 


3 


for) 
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Stability region for B)>0, B,>0 , B,+B,=1 


scaled by N 


-6 -4 -2 0 2 
scaled by N 


Figure 3: The stability region of the new RK method with Bo > 0,61 > 0,89 + 61 = 1 


In Figures 1, 2, 3, the absolute stability regions of the given RK method 
(11) have been demonstrated for various values of 89 and 6,. We observe 
that in the case 89 > 0, 8; > 0,89 + 61 = 1 the new RK method is A-stable 
and in the cases 6) = 0,6, = 1 and (o = 1,8, = 0 the absolute stability 
regions, in terms of h, tend to the same regions of the forward and backward 
Euler methods, respectively. 


4 Linear multistep methods 


In this section, we present the general theory for an arbitrary linear multistep 
method. We demonstrated the main idea by using the one-step and multi- 
step starting procedures. In both cases, we obtain the corresponding Runge- 
Kutta scheme and the LTE of this method provide an estimation of the global 
truncation error of the given LMM. 


4.1 A single step method as starting procedure 


Now, we represent a general linear multistep method (LMM) in the form of a 
Runge-Kutta scheme. For simplicity, we use the trapezoidal rule as starting 
procedure of the LMM. We show the starting values {yn+,; a as a linear 


combination of the y,, and { peers Therefore, starting with 


108 J. Farzi 


h 
Ynt1 = Yn + 5 (in + Fn+1)s 


we find 


h 
Yn+j = Yn t 5 (fn + Qfngi te: +2fngj-1t fri), J =1,2,... (19) 


therefore, we have 


k-1 k-1 


{( aj)fn t+ (ar+25— 05) fri 


k k-1 
s; AjYn+j = Yntk + Oo 5) Yn + 
j=0 j=0 j=l jee 


+++++ (Qp—2 + 20%-1) fn4e—2 + Qr—-1fnte—1)}- (20) 


Substituting (20) into the LMM (1), we obtain 


k-1 k-1 
Yntk = — y OjYn+5 + hy- Bs Fai 
j=0 j=0 


k-1 1 eat 1 k-1 
= =; 045 )Yn + AL (Bo — 3 S05) fa aa ie go = So aj) fn+1 
j=0 j=l j=2 


1 
tot (Boz = 5-2 — On-1) Fn te—2 


+(Be—-1 — Sail ones + Befn+e)}, 


2 
or, 
k+1 
Untk =Yn th >) bifnsi-t, (21) 
i=1 
where, 
ye 
br =o 5 Ds Qj, 
j=l 
1 k-1 
bi = Bia — 501-1 — Day, 2<i<k-l 
j=i 


1 
be = Be-1 — gek-1 
be+i = Gr. 


It is straightforward to show that 
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On defining 


. h 
Kiai = f(€n46) yng) = flan +h, Yn 5 (Fn EZ Feist ef, iia hone)» 


the scheme (21) can be represented in the form of the following RK method 


k+1 


Ym+1 = Ym + Hy vK:, 
i=1 


where, 


We will change the subscript n to m in order to show that when LMM runs 
from Yn to Yn+x~ the RK scheme runs just a single step from Ym = Yn to 
Ym+1 = Yntk- Therefore, we set 


tm = Lnytm41 = In + A, 


aij= 4p 255G<i-1 


a 
Gq = k” 
b; 
= — 
7 k’ 
H=kh. 


In this case, the first order condition holds 


k+1 


oh 
w=1 


but, the second order condition no longer holds 
i=1 ier ae 


In k > & this condition turns out to be the exact order condition. Again, 
the principle local truncation error (PLTE), where p = 1 reads 
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H? 
PLTE = > 


& 2Ol- 1091 = — 5 h(x — 20)y) (20) 


and therefore, we have 
1 
LTE = —5h(ae — ao)y™ (ao) + O(h?), 


which is the global error in k times application of a k-step LMM with Trape- 
zoidal rule as starting procedure. 


4.2 A Runge-Kutta method as starting procedure 


Now, we have considered the general explicit Runge-Kutta method (5)-(6) 
as the starting procedure for the linear multistep method and we find the 
Runge-Kutta representation of the given LMM (1). Based on this starting 
procedure, we have obtained approximate values for yn+4;,j =1,...,k—las 
follows 


s 
Ynt+j = Yn +jh>- bik\?), 
i=1 
Ko = f(tn + jch, yn + jh S~auk”), i=1,...,s+1, 
1=1 
where ¢s41 = 1, @i,s41 = 0,4 =1,2,...,8+1,@541,; = 6;,7 =1,2,...,8. For 
an implicit method corresponding to the yn+%, we define 


s+1 
Kh) = f(an + kes4ith, yn + khS- asi sk). 
l=1 


By inserting these approximations into the LMM (1) we obtain 


k-1 k 
Yn+k = — » A5Yngj + oe Bi fn+s 


k-1 s 


Dah — AS SS ajbik? en doant ee 


g=l1 t=1 


I 


there is s(k —1)+ 2 different KY ) in the above representation, however we do 
not distinguish them and consider (s+1)(k—1)+1 moments. The advantage 
of ignoring the similarity in the moments is that the resulted Butcher array 
is simpler to work and it is convenient to prove the theorems. Let, 
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tm = %n,H =kh,tmy41 =tm + H = Unik 
and 
Cit (G—-1(st1) = La, a=1,...,8+1,7 =1,...,k-1, 
a 1 
C(s+1)(k—-1)4+1 = rm 
and 
ae a — = (arb Bo): 
py) = —Fjaybs, é=1,...,8, 7=1,...,k-1, 
p= Bind =1,...,k, 
the vector b consist of these values: 
B= [BM ,b2), 2. bO-D, o™ ]. 
Therefore, we find the new Runge-Kutta scheme 
Ym+1 = Ym + HY biki, 
i=1 
i) j ee Oe. op 
ke? = f(tm + FeH, Ym + p> ak; ), t=1,...,8 41, 


i=1 
where, § = (s+1)(k-—1)+1 and 


= k 
ong RO ROO, 


The corresponding Butcher array is given in Table 1. where 0 is a (s +1) x 1 
zero vector. Introducing D as a (k — 1) x (& — 1) diagonal matrix 


The more compact form of the above scheme is resulted. 
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Table 1: Butcher array of the Runge-Kutta representation of LMM (1) scheme 


A 
Le t[7o| 
2 ah 
= cae tA 
bed [$0 
k pr 
pr 


Table 2: Butcher array of the Runge-Kutta representation of LMM (1) scheme 


nef 


pr 


pr 


Now, to verify the order of the new Runge-Kutta scheme suppose that the 
order of LMM (1) and Runge-Kutta scheme (5)-(4) are p and #, respectively. 
Then, we can prove the following theorem. 


Theorem 2. If the Runge-Kutta scheme as starting procedure has order 
p and the order of linear multistep method (1) is p, then the corresponding 
Runge-Kutta method with Butcher array in Table 2 is of order min{p, p}. 


Proof. The order condition corresponding to the tree 
t NY 2 


with m leaves is 


= 1 
‘ii — cand 
2, bit aa m=0,1,...,p 


for any m < min{p, p} we prove that 


*. 1 
De aa m =0,1,...,min{p, p}. 
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According to the definition of @ and b, we have 


5 k-1 s ; k 
7 =m 1 . J m 1 J m 
5 bey" = — 5 bee) > bi(Fci)™ + if > Bi(Z) 
i=1 j=l i=l j=l 
k-1 k 
1 1 1 
= m+1. “M2. 
~  m+1km41 23 a o Emel di 8; 


k k 
1 1 
_ m+1 -m+1. M2. 
~ ee (Oa) + 
j=0 j=0 


1 
m+ 


note that in the above relations we have used the order conditions for starting 
Runge-Kutta method as well as the order conditions for linear multistep 


methods: A 


k 
a ita, = FB in = 204 


j=0 j=0 


The rest of the proof is closely related to the block structure of the Butcher 
array in Table 1. The extracted elements of the Butcher array is demon- 
strated in the Table 3. 


Table 3: Elements of Butcher array of Table 1 
A 
l l 


Ee | |p? 
; 1=1,2,...,k-1 
yp) 
1 br 
= = GPx 


These partial elements will help us to exactly find the effect of them in 
multiple sums in order conditions. 

_The last condition to complete the proof of the third order conditions is 
>. jel bid c;. The role of each element in summation, separately, is 


oa *eubi)(4-4%5)(3) (A) Derg) i=1,2,...,k-1 


ij=l i=1 


Pal be) 
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summing up these terms we obtain 


k=1 s s 
De aves = Li (Me peda Geant gar + GH) NGG) 
+o bie) 
1 1X 
= GR (k8 — $0 Bay) + 553 SOPH 
7 1 I=1 1=0 
6 


similarly we can prove the higher order conditions. 


4.3 Local truncation error of the new Runge-Kutta 
method 


We have shown that the order of the new Runge-Kutta method with Butcher 
array in Table 2 is p* = min{p,p}. Ignoring the effect of the roundoff errors, 
we can consider LTE of this method as the global error of the given linear 
multistep method in evaluation of yn+,. The LTE of this scheme now reads 


LTE = Gen S> aft)[L- 1@QvWIFO + O(n? *). (22) 


5 Numerical illustrations 


Example 1. As an example we consider the Heun’s third order 3-stage 
formula 


WIN WI © 


dj 
3 

2; 
05 
T 
4 


as starting procedure for the third order convergent linear multistep method 


1 1 3 h 
Yn+3 + qunte nti gun = 8 [19 fr+e ee 5 fn], (23) 


the corresponding Runge-Kutta method is 
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0 

ifa 

9| 9 

2 2 

9/9 9 

1/1 

3|i3 9 ia 

010 000 iui 

2 2 

2020:-0: 0. 2 

4 4 

21070200: 0° 4 

2 2 6 

2).0°0030° 2 

1 si 0 i 0 Hi 0 si 5] 0 
xi 0 i 0 Hi 0 3 51 0 


This method is of order p* = min{3,3} = 3. Substituting the above data in 
LTE (7) we find that 


4 
ec ‘ S> a(t)[L — y(t) dO)|F + O(h°) 
r(t)=4 
4 
=F (Bu - Bera Bots) +00), 


where all functions are evaluated at x = x, and y = yy. Note that this 
formulation maintains the actual order of both schemes, while the error term 
is only exact for {1,2}. 


To numerical illustrations, we consider (2) with the following data [9] 


f(x,y) = [v, v(v = 1)/uJ’, LE (0, 1], (25) 
where, 
y= [u,v]*, y(0) = [1/2, —3]". 


In this test we take N = 51 with h = 1/50 for 3-step method (23) and 
H = Nh = 3/50 for the corresponding Runge-Kutta scheme (24). 


Figure 4 illustrates the global error (accumulation error) of (23) for test 
problem (25). The estimation of GTE of (23) is shown in Figure 5. As 
we have proven the LTE of (24) is the GTE of (23). However, to find the 
true LTE we make localizing assumptions in implementation of (24), i-e., in 
evaluation of y,4+1 we assume that y, = y(x,). The comparison of the third 
portions of Figure 4 and Figure 5 justify the efficiency of the given estimation. 
The negligible difference in the error is due to roundoff errors. 
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lull, 
: 


Figure 4: The numerical (red circles) and exact solutions (solid line) of (25) with (23), 
and the GTE of the method 


u 
Oo 
wo 

Vv 

1 
4 
oa 

lull, 


0.1 : a : 0 
0 : 0 


Figure 5: The numerical (red circles) and exact solutions (solid line) of (25) with (24), 
and the LTE of the RK method 


Example 2. In this example, we consider a four-step, fourth-order total 
variation bounded (TVB(4,4)) linear multistep scheme (1) with the data 
given in Table 4 [11]. 
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The coefficients of the four-step, 
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fourth-order linear multistep scheme 


ay 


Bi 


Table 4: 

(TVB(4,4)) 
1=0 
t=1 
1=2 


t=3 


(= 


—0.345464734400857 
1.494730011212510 
—2.777506277494861 
2.628241000683208 
1 


—0.620278703629274 


2.229909318681302 


—3.052866947601049 


1.618795874276609 
0 


We utilize the fourth order classic 


procedure. 


It turns out that the corresponding 
in Butcher array (27), where the el 


| NIP NIF © 
alrFl OC Onik © 
we! © NIF 


wle| rR 


The elements of vector b in Runge-Kutta scheme (27) 


bj 


bits 


Table 5: 
i=1 
i=2 
1=3 
i=A4 
1=5 
1=6 
ee 
i=8 


—0.092789258773464 
0.124560834267709 
0.124560834267709 
0.062280417133855 
0.557477329670325 

—0.231458856457905 

—0.462917712915810 

—0.462917712915810 


—0.231458856457905 
—0.763216736900262 
0.328530125085401 
0.657060250170802 
0.657060250170802 
0.328530125085401 
0.404698968569152 
0 


Runge-Kutta method as a starting 


(26) 


Runge-Kutta scheme takes the form given 
ements of b is given in the Table 5. 
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Ye OC Oar © 


Slr O or 


Blk OC Oar Oo 
alk © le 

lH wie 

m= 

RIF 


a OO aw © 
Bl CO lw 
AIR loo 


hr Blew iw colw colwW OO NIP NIP Al Ale © AIR AIR GIR OIF © 
Cole 


b 
b 


The test problem is the same as the previous example. For numerical il- 
lustrations we take N = 120, thus we have h = 1/120 and accordingly 
H = 4/120 = 1/30. The numerical results including the GTE of TVB(4,4) 
scheme and LTE of (27) are presented in Figure 6 and Figure 7, respectively. 
By comparison, we find the excellent estimation of global truncation error of 
TVB(4,4) based on new Runge-Kutta method. The maximum error in this 
estimation is 2.67& — 05. 


6 conclusion 


In this paper, we have developed an estimation for the global truncation error 
of a linear multistep method. The global error analysis is more complicated in 
comparison with the local error analysis. The key idea is the representation 
of several steps of the LMM as a single step of a corresponding Runge-Kutta 
method. Therefore, the analysis of global error of a LMM accomplished by 
estimating the local truncation error the corresponding new Runge-Kutta 
method. We have demonstrated the theoretical aspects for some important 
class of linear multistep methods with total variation bounded (TVB) prop- 
erty, which is a crucial property in selecting an appropriate time marching 
method for solving nonlinear conservation laws [11]. 
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Figure 6: The numerical (circles) and exact solutions (solid line) of (25) with TVB(4,4), 
and the GTE of the method 


x10 
0.5¢ 2.5 
0.45 
2 
0.4} 
0.35 5 15 
_N 
5 4 > -1.5 = 
1 
—-2H 
0.5 
-2.5 
—36- 0 
0 0.5 1 0 0.5 1 
xX xX xX 


Figure 7: The numerical (circles) and exact solutions (solid line) of (25) with (27), and 
the LTE of the RK method 
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